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Abstract 

Particle methods are widely used because they can provide accurate 
descriptions of evolving measures. Recently it has become clear that 
by stepping outside the Monte-Carlo paradigm these methods can be of 
higher order with effective and transparent error bounds. A weakness of 
particle methods (particularly in the higher order case) is the tendency for 
the number of particles to explode if the process is iterated and accuracy 
preserved. In this paper we identify a new approach that allows dynamic 
recombination in such methods and retains the high order accuracy by 
simplifying the support of the intermediate measures used in the itera- 
tion. We describe an algorithm that can be used to simplify the support 
of a discrete measure and give an application to the cubature on Wiener 
space method developed by Lyons, Victoir [12) . 



1 Introduction 

In pricing and hedging financial derivatives as well as in assessing the risk inher- 
ent in complex systems we often have to find approximations to expectations of 
functionals of solutions to stochastic differential equations (SDE). We consider 
a Stratonovich stochastic differential equation 

d 



* Research supported by a grant of the Leverhulme trust {F/08772/E). 
t Original version: June 7, 2007 



1 



defined by a family of smooth vector fields Vi and driven by Brownian motion. 
It is well known that computing Px-tf ■= ^'l/lCr-t.^)) corresponds to solv- 
ing a parabolic partial differential equation (PDE). High dimension and hypo- 
ellipticity are common obstacles that arise when one calculates these quantities 
numerically. When facing these obstacles some classical computational methods 
become unstable and/or intractable. 

There are many settings where one is interested in tracking the evolution of 
a measure over time in an effective numerical fashion. One example is the nu- 
merical approximation to the solution of a linear parabolic PDE. In this case one 
tracks the evolution of the heat kernel measure associated to the PDE. Another 
example is the filtering problem where one wishes to approximate the unnor- 
malised conditional distribution of the signal, which is governed by a stochastic 
partial differential equation known as the Zakai equation. 

An evolving measure can be viewed as a path in the space of measures. Thus, 
even if the underlying state space is finite dimensional, we potentially face an 
infinite dimensional problem. Particle approximations can in many cases pro- 
vide good descriptions of evolving measures, see for example the survey articles 
[21 13] • Higher order methods may allow us to take far fewer time steps than 
classical methods in the approximations. An example of a higher order particle 
method may be found in Kusuoka [5]. Although effective in practise (compare 
Ninomiya [M] and Ninomiya, Victoir [13]), these methods have the drawback 
that the number of particles can explode exponentially if the process is iterated 
and accuracy preserved, see for example Lyons, Victoir |12] . 

Sometimes the essential properties of a probability measure we care about can 
accurately be described and captured by the expectations of a finite set of test 
functions. If we can find such a family of test functions we can replace the 
original measure with a simpler measure with smaller support that integrates 
all test functions correctly and hence still has the right properties, provided of 
course the number of test functions is small compared to the cardinality of the 
support of the original measure. We will also insist that the reduced measure 
jl has supp{jl) C supp{ii). This condition ensures that feasibility constraints 
imposed on the measure ^ will also be satisfied by fl. For a finite Borel measure 
on a pohsh space and a set of integrable functions {pi, . . . ,p„} we can show 
that such a reduced measure jl always exists with card{supp{jl)) < n + I . 

In this paper we present a simple algorithm that can be used to compute re- 
duced measures for discrete measures ^. The runtime is polynomial in the size 
of the support of the measure ^. The algorithm relies on the observation that 
if P is the i?" valued random variable P{x) :— {pi{x), . . . ,Pn{x)) and fip the 
law of P under the measure /i, then finding a reduced measure fl is equivalent 
to finding flp a discrete measure on i?" with card(supp(/ip)) — n + 1 and the 
same centre of mass (CoM) as /ip. 
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We describe an application to the Kusuoka- Lyons- Victoir (or KLV/cubature on 
Wiener space) method developed by Lyons, Victoir in [12], following Kusuoka 
[S]. It provides higher order approximations to E{f{^j,^)), if the test func- 
tion / is Lipschitz and the vector fields satisfy Kusuoka's UFG condition (see 
[7]) which is weaker than the usual Hormander condition. The expectation 
^ifi^Tx)) niight be viewed as an infinite dimensional integral against Wiener 
measure. The authors construct discrete cubature measures Qt = X)r=i ^i^i^i.T 
supported on continuous paths of bounded variation that approximate Wiener 
measure in the sense that they integrate iterated integrals up to a fixed degree 
correctly. The expectation of a Wiener functional f{£,Tx) against the discrete 
cubature measure may be obtained by computing the endpoints of the solu- 
tion of the SDE along the paths in the support of Qt- Thus the KLV method 
might be viewed as a discrete Markov kernel taking discrete measures on to 
discrete measures on . More explicitly we have 

n 
i=l 

and 

^QtH^T/x) = EKLV{S^,T)f- 

The bound on the error when replacing the Wiener measure with a cubature 
measure is given in terms of higher order derivatives of /, so in general will not 
be small as / is only assumed to be Lipschitz. The results in Kusuoka, Stroock 
[B] and Kusuoka 7^ show that Ptf will be smooth, at least in the direction 
of the vector fields Vi. This is resolved by applying the method iteratively 
over a partition of the time interval [0,T]. The operator corresponding to the 
iterated application of the KLV method is Markov and hence, the error of the 
approximation of Prf on the global time interval [0,T] is the sum of the error 
of the approximations over the subintervals of the partition. So considering 
an uneven partition of the global time interval [0,r] with time steps getting 
smaller towards the end, we can iteratively apply the cubature method over the 
subintervals and reduce the error in the approximation to any accuracy. If m is 
the degree of the cubature formula we can find a partition such that the error 
in the weak approximation is uniformly bounded by 

Cfc-(™-l)/2||V/|U, 

where k is the number of time steps in the partition and C a constant indepen- 
dent of k and /. 

The iterated KLV method might be viewed as a particle system on R^ where 
the particles branch in an n-ary tree. Hence, the number of ODEs to solve grows 
exponentially in the number of iterations. In this paper we add recombination 
to the KLV method. After each application of the KLV operation we replace the 
intermediate measures by reduced measures. The property of the KLV measure 
we are targeting is to integrate Ptf, the heat kernel applied to /, correctly. We 
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have identified a finite set of test functions that ensures that the bound on the 
overall error of the approximation of Pt/ is only increased by a constant factor 
and hence the modified method has the same convergence properties. Moreover 
we can show that under the Hormander condition for bounded vector fields the 
number of test functions required grows polynomially in the number of itera- 
tions. 

We believe that the combination of the two ideas - higher order particle meth- 
ods to describe the evolution of a measure on the one hand and simplifying 
the support of the measures used in the description, by characterising essential 
properties of a measure using the expectations of a finite set of test functions on 
the other hand - have more general applications than investigated so far. Appli- 
cations to the stochastic filtering problem appear to be particularly promising, 
see Litterer, Lyons [5], [THj for an outline. 

2 A reduction algorithm for the support of a 
discrete measure 

Let us start the precise description of the reduction problem. The notation 
in this section is independent of the notation used in the description of the 
cubature method in the following sections. Consider a finite set of test functions 
Pn = {pi, . ■ ■ ,Pn} on (ri, /x), a measure space with /i a finite discrete measure 



with large support. By this we mean that h is at least of order n^. In the fol- 
lowing we assume that /x is a probability measure, i.e. the weights add up to one. 

Definition 1 We will call a discrete probability measure p, a reduced measure 
with respect to /i and P„ if it satisfies the following three conditions: 

1. supp{li) C supp^ji) 

2. For all p & Pn 



3. card{supp{jl)) < n -I- 1. 

The first condition is more important than it looks as it ensures that feasi- 
bility constraints imposed on samples drawn from fi will also be satisfied by fl. 
We wish to construct effective algorithms to compute the reduced measure. 





i=l 
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Let P be the i?" valued random variable P := (pi, . . . ,p„) defined on (fi, /i). 
Then the law jip of P is the discrete measure on i?" 

n 

/ip = ^Ai5:j, = (pi(zi), . . . e i?". (1) 

1=1 

The centre of mass (CoM) for the measure /ip is given by 

n 

CoM(Aip) = ^A,a:,. (2) 

1=1 

To find a reduced measure we articulate an equivalent problem in terms of /Zp. 
The problem becomes to find a subset Xi,^ of the points Xi and positive weights 
Xif. to produce a new probability measure jlp such that CoM{Jlp) = CoM{pp). 
A reduced measure jl is then easily obtained from jlp by taking 

with Zij^ e supp{fj,) satisfying P{zi^) = Xi^. 

Note that given any subset Xi^, there exist suitable weights A^^. if and only 
if CoM{^p) is contained in the convex hull of these points. Caratheodory's 
theorem implies that in principle one can always find jlp with support having 
cardinality at most n + 1 and the algorithm explained below provides a con- 
structive proof to that. 

By considering Xi — CoM(fj,p) in place of the Xi we may assume without loss of 
generality that CoM{fj,p) is at the origin. We may also assume that the Xi are 
all distinct, as we can otherwise eliminate points Xi from the original measure 
H by sorting and combining them. 

A first algorithm (Algorithm 1), communicated to us by Victoir [T5], sequen- 
tially eliminates particles from the support of the measure. It is well known 
and has for example been used in constructive proofs of Tchakaloff's theorem 
(Davis 0). 

Given any n + 2 points, the system given by 

n+2 

UiXki = (3) 

1=1 

n+2 
i=l 

is a linear system with n + 2 variables, but only n + 1 constraints. Therefore it 
has a non-trivial solution, which may for example be determined using Gaussian 
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elimination. Thus we may either add 



mm 

Ui<0 




to ([2]) or subtract 



mm 

Ui>0 




from ^ leaving all weights in the result non-negative and their overall sum 
unchanged. In either case, by construction, the coefficient of some Xj vanishes. 
We now have obtained a new probability measure with the same centre of mass 
and at least one point less in the support. Applying the procedure iteratively 
until there are only n -\- 1 points left we obtain a reduced measure. Clearly the 
method requires no more than n iterations of the above procedure. 

Remark 2 If h is the dimension of the lowest dimensional (affine) subspace 
of i?" containing the set {{pi{y), . . . ,Pn{y)) \ V € supplp)}, we can continue to 
apply the elimination procedure described in Algorithm 1 until card{supp{fl)) < 



For improving the order of the overall algorithm we now look at suitable 
linear combinations instead of points. 

To describe the algorithm we define an abstract Procedure A that takes a dis- 
crete probability measure v with 2(n -I- 1) particles in its support and returns 
another discrete probability measure v with (n -\- 1) particles in its support 
satisfying CoAI{v) ~ CoAI{v) and supp{i>) C supp{v). Procedure A may for 
example be realised by n -I- 1 applications of the reduction procedure of Algo- 
rithm 1. 

Main reduction algorithm (Algorithm 2): 

1. Partition the support of /ip — X]r=i ^i^xi into 2{n + 1) sets of as near equal 
size as possible. Let these sets be denoted by Ij, 1 < j < 2{n + 1). 

2. Compute the probability measure v = X]i=i^^'' ^i^xi where 



h + 1. 



and lyj = fip{Ij) = Y.^■.x,<^l, 

3. Apply Procedure A to compute a measure v — "Yli 



.n+l 

'i=i 



Vi^&i^^ with CoM(zy)=CoM(i>) 



4. Repeat 1.-3. with 



i^p = Y^ ""^.T^^-^ 



for /ip until n -f- 1 particles are left in the support of /ip. 
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Proposition 3 Given ^ and Pn the algorithm described above requires [lg(n/n)] 
iterations of Procedure A to compute a reduced measure. 

Proof. We might interpret the points Xj as the respective centre of masses of 
the individual subsets Ij. 

It is clear that fi'p has positive weights and support contained in the support of 
/ip. Hence, we only need to show that CoM{pi'p) = CoM{pip). 

We have 

3 = 1 XkGlij j=l 
2{n+l) 2(ri+l) 

= CoM{i)) = CoM(v) = Vjij = Y ""^Y ~ = CoM{fip). 

As n < n2 r's("/n)l ^lay assume without loss of generality that n = n2 ^ 
It is obvious that each iteration halves the number of particles in the support 
of fj,p and we require exactly [lg(n/n)] iterations. ■ 

Corollary 4 Using the main reduction algorithm we can compute a reduced 
measure with respect to fi and Pa in 

0{nn + n \og{n/n)C{n + 2, n + 1)) 

steps where C(n + 2,n + 1) represents the number of steps required to solve a 
system of linear equations with n + 2 variables and n + l constraints. 

Proof. To compute the intermediate measures v we need to calculate n- 
dimensional linear combinations. The number of steps required for these addi- 
tions is bounded above by the series 



oo 

n 

i=0 



^n2-* = 2r 



The procedure A may be realised by n+1 applications of the reduction procedure 
used in Algorithm 1 described above. ■ 

Remark 5 Note that the linear systems of equations we need to solve in the 
algorithm are singular. Hence, for a practical implementation we have used a 
method based on the singular value decomposition (SVD) to avoid numerical 
instability. 



dll with an implementation of a version of the algorithm and a Visual 
Studio project with a simple example for its use can currently be found at 
[http: / /www. maths. ox. ac.uk/~tlyons/Recombination/reduce-dist_01^_ paper. zip 
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If the support of the measure /j, we wish to target is particularly large or 
possibly even infiiiite we can consider a different approach. If we can find a 
subset of points that with a reasonably high probability contains the CoM in 
its convex hull, we may use linear programming to check if a given set of points 
contains the CoM in its convex hull and reconstruct the weights. The results 
in Wendel [TB] imply for example, that a collection of k uniform IID random 
variables on the unit sphere in contains the origin with probability 



j=o 

In particular this yields Pn,2N — 1/2. 



k-l 

i 



3 Outline of the cubature algorithm 

We describe the cubature method developed by Lyons and Victoir [12] . Through- 
out the paper C is a constant that may change from line to line, specific con- 
stants however will be indexed C1C2, ■ ■ ■ ■ Let C^{R^ , R^) denote the smooth 
bounded R^ valued functions whose derivatives of any order are bounded. Then 
V^ e C^{R^,R^), <i<d may be regarded as vector fields on i?^. We de- 
fine a partial differential operator L = Vq + \ [Vi -t- . . . V]i) and consider the 
following parabolic partial differential equation (PDE) 

du 

-Trr{t,x) = -Lu{t,x), 

dt (4) 

u{T,x) = f{x) 

for a given Lipschitz function /. The aim is to find an approximation of 
u{0,x) for a given x. Consider the probability space (Cq ([0, T], i?''), J", P), 
where Cq ([0, T], i?'') is the space of i?'^ valued continuous functions starting 
at 0, T its usual Borel tr-field and P the Wiener measure. Define the coor- 
dinate mapping process BKui) — uj^{t) for t £ [0,r], a; e O. Under Wiener 
measure, B = {B} , . . . , Bf) is a Brownian motion starting at zero. Furthermore 
let B^{t) = t. Let ^, t G [0,7^], x € R^ be a version of the solution of the 
Stratonovich stochastic differential equation (SDE) 

d 

di,^, = Y,V.{i,^,)odB\A^,. = x (5) 

4=0 

that coincides with the pathwise solution on continuous paths of bounded vari- 
ation. In this case classical theory tells us that u{t,x) = E{f{£,T-t x)) is the 
solution to (U]). 

We define the Ito functional $t,x : C^{[Q,T], R"^) -> R^ by 

<i>T,.(w)=eT,.(^)- (6) 
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Denote by Rm [Xi , . . . , Xd] the space of polynomial^ in d variables having de- 
gree less or equal to m. Let be a positive Borcl measure on W^. A discrete 
measure jl 

n 

with xi, . . . ,Xn contained in supp{fi) satisfies a cubature formula of degree m if 
and only if for all polynomials P £ R„i [Xi , . . . , Xj] , 

f P{x)fl{dx)= f P{x)ll{dx) \P{Xi). 
JR-^ jR-i ~{ 

It is well know that if all moments of /i up to degree m exist we can can 
always find such a measure with 

card(supp(/i)) < dim(R,„[Xi, . . . , Xd]) + 1, 

see e.g. Bayer, Teichmann [1 . More generally we have the following lemma, 
which we state without proof. 

Lemma 6 Let be a polish space, T its Borel sets and fj, a Borel probabil- 
ity measure on (fi, !F). Let /i, . . . , /„ be a finite sequence of real valued Borel 
measurable functions on the probability space {Q,,J-,^) with E{\fi\) < oo for 
1 < i < n. Moreover suppose that D is a Borel set with fJ-{D) — 1. Then there 
exist points wi, . . . , Wn+i € D and a discrete measure 

n+l 

i=\ 

such that 

E^{fO=E~^[f,) 

for \ < i < n. 

In other words [i admits a reduced measure fip with respect to any finite 
set P of integrable functions. In connection with the use of the Taylor formula 
a cubature measure provides an effective tool for integration over finite dimen- 
sional spaces. 

One can formulate an analogous condition to identify cubature measures on 
Wiener space. Here the role of polynomials is taken by iterated integrals of the 
form 

/ odB\\---odB\l. 

JO<ti<-<tk<T 

We identify this iterated integral by the multi-index (ii, . . . , ife). 

^Any finite dimensional space of integrable and continuous functions could be used to 
define cubature. This extension can be helpful 
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Define the set of all multi-indices A by 

oo 

A= |J{0,...,4'= 

and let a — (ai, . . . , afe) € A be a multi- index. Furthermore we define a degree 
on a multi-index a by ||q:|| ^ k + card{j : Uj — 0) and let 

A{j) = {a e A : \\a\\ <j}. 

Moreover define Ai by Ai = A\{0, (0)} and let Ai{j) ^ {a e Ai : \\a\\ < j}. 
It follows from the scaling property of Brownian motion that 

odSfi ■■■odB?"'. 

tl Cfc 

/0<ti<-<tfc<T 

equals, in law, 

rii"ii/2 [ odB-^l ■ ■ ■ o dB^^ . (7) 

"'0<ti<-<tfc<l 

Definition 7 Fix a finite set of multi-indices A Q A. We say that a discrete 
measure Qt assigning positive weights Ai, . . . , A„ to paths 

wi,...,w„eCoV([o,r],i?'^) 

is a cubature measure, if for all (ii, . . . , i^) € A, 



\JO<ti<-<tk<T / ^^-^ JO<ti<-<tk<T 



duHtk) 



where the expectation is taken under Wiener measure. If A = A{m) we say that 

n 

is cubature measure of degree m. 

In |12| . the authors show that one can always find a cubature measure sup- 
ported on at most card(yl) continuous paths of bounded variation. More im- 
portantly they give an explicit construction of a degree 5 cubature formula with 
O(d^) paths in its support. 

Suppose paths wi, . . . ,a;„ and weights Ai define a cubature measure for T — I. 
It follows immediately from ([7]) that the measure supported on paths ujt i given 

by 

ujir,,^VTu;lit/nj = l,...,d (8) 
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and unchanged weights defines a cubature measure for general T. From now 
on suppose that the measure Q := Qi is a cubature measure of degree m. 



The following proposition, taken from jI2j . is the key step in estimating the 
error Et when one approximates the expectation of f{^Tx) under the Wiener 
measure by its expectation against Q. 



Proposition 8 



Et '■— sup 



m+2 



/2 



sup 



IK 



j=,n+l (ai,...,a,)GA(i)\A(i-l) 

where C is a constant that only depends on d, m and Qi. 

In general the right hand side of the inequality in Theorem [8] is not sufficient 
to directly obtain a good error bound for the approximation of the expectation, 
in particular if / is only assumed to be Lipschitz the estimate appears useless. 
So, instead of approximating 

Prfix) Eifi^^J) 

in one step, one considers a partition V of the interval [0,T] 

to = < ii < . . . < tfc = T, 

with Sj = tj —tj-i and solves the problem over each of the smaller subintervals 
by applying the cubature method recursively. If r and r' are two path segments 
we denote their concatenation by t ®t' . For the approximation we consider all 
possible concatenations of cubature paths over the subintervals, i.e. all paths of 
the form w.,, ® . . .®u 



, . We define a corresponding probability measure v 



by 



Ail • • • ^ik^LO,^,i 



4l,...,lfc = l 

The following theorem taken from Lyons, Victoir |12j is the main error estimate 
for the iterated cubature method, which we in the following also refer to as the 
Kusuoka-Lyons- Victoir (KLV) method. 

Theorem 9 The total error Ex> for the approximation 

Ej) 



sup |PT/-i?.(/(eT,J)| 



sup 



'^Sk,^k)) 



11 = 1 
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is bounded by 

(m+l fc-1 (j+l)/2 \ 
j—m i—1 j 

where C\ (T) is a constant independent of f and k, the number of time steps in 
the partition of the time interval [0, T] . 

To compute the expectation with respect to the measure v exactly requires 
one to solve 

n^+i - 1 
n — 1 

inhomogeneous ODEs - each corresponding to a path ijJsi.ii ^ ■ ■ ■^'^Sk,ik ~ where 
n denotes the number of paths in the support of the cubature measure Q and 
k the number of subintervals in the partition. Hence, the number of ODEs to 
solve grows exponentially in the number of iterations. 

Following Kusuoka [7] we define for multi-indices a = (ai, . . . , a^), /? = 
{Pi, . . . , f3i) & A a multiplication by 

a*l3= {ai,...,ak,f3i,...,Pi). 

We inductively define a family of vector fields indexed by A by taking 

V[0] = 0, =V^, 0<i<d 

V[a„] = [V[a], K], 0<i<d,ae A. 

The main ingredients used when obtaining the bound (jH]) are Proposition [5] and 
the following regularity result due to Kusuoka and Stroock 6 and Kusuoka [7] , 
which says that even if / is not smooth Psf is smooth in the directions of the 
vector fields Vi. Let / be Lipschitz and ai, . . . ,ak £ Ai then for all t E (0, 1] 

11^..] . . . ^„.]PJI|oo < ,(||a.|| + ...+ |K||)/2 llWHoo (10) 

provided the vector fields satisfies the UFG condition defined below. 
Following Kusuoka we introduce a condition on the vector fields 

Definition 10 The family of vector fields Vi, i — 0,...,d is said to satisfy 
the condition (UFG) if the Lie algebra generated by it is finitely generated as a 
left module, i.e. there exists a positive k and Ua.fs G Cj^ satisfying for all 
a G Ai, 

V[a] = ""./3^[/31- (11) 
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The bounds for the error of the KLV method derived in Theorem [S] (see 
Lyons, Victoir [12] for details) assume that the system of vector fields Vi, 
i — 0, . . .d satisfies the UFG condition. 



Definition 11 We define the (formal) degree of a vector field Vj^jja G A de- 
noted by da to be the minimal integer k such that may be written as 

/3eAi(fe) 

with Ua^p € C^. 

Note that for a G Ai we always have da < It was pointed out in Crisan, 
Ghazali [5] that the analysis in Lyons, Victoir [T^] for the bound in © requires 
Vq to have formal degree at most 2. If the formal degree of Vq is greater the 
bound in changes and all bounds in the paper will change accordingly. For 
sake of simplicity we will in the following assume that Vq has formal degree 2. 
The bounds can be improved in an obvious way if the degree is 1 or 0. For a 
generalised error estimate based on Kusuoka's ideas ([S]) that does not require 
this additional condition see Litterer [TT| . 

A trivial generalisation of Corollary 18 in Crisan, Ghazali [3] allows us to 
state a version of the Kusuoka and Stroock estimate in terms of the formal 
degree of a vector field. Let / be as above and ai,...,afc e A then for all 
t e (0, 1] 

\\V[a.] ■ ■ • V[a,]Ptf Woo < ^(,,^+...+,,^)/J |V/Hoo. (12) 

For the remainder of the paper when we consider recombination we are going 
to assume the following uniform Hormander condition. 

Definition 12 We say that a collection of smooth vector fields Vi, i — 0, . . . ,d 

satisfies the uniform Hormander condition (UH) if there is an integer p such 
that 

infJ (V[„](a;),0';a;,eei?^,|^| = li :=A/>0. 

Note that the uniform Hormander condition implies the UFG condition. 
Under this stronger assumption it is straightforward to show in addition that 
Ptf is a smooth function on with explicit bounds on its derivatives. We 
outline an argument below that follows Kusuoka [7 and gives bounds on the 
regularity of Ptf, which we will use in the following section when we apply 
recombination to the cubature method. 
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Following Kusuoka fT" let F{x) G C^{R^:, (E) R^) be given by 
F{x)= J2 ^a](^)®^H(a;), xeR^, 

and Aq : R^ — >■ [0, oo) be the continuous function 

Xo{x)=mf{{F{x)y,y);yeR'' ,\y\ = l}, x e R^" . 

Note that 

{F{x)y,y)= J2 

aeAi{p) 

and hence under the Horniander condition (UH) we have Ao(a;) > ^/ > 
for all X £ i?^. As in Kusuoka iT] let = and Ua^i : i?^ R, 

a e Ai{p),i — 1 , . . . , iV, be given by 

ac^Ax) = (e„ f^(a:;)-il^„] (a;)), x € R^ (13) 

and observe that 

d 

— = (^c.,^V[o.]. (14) 

The following lemma may be found in Kusuoka T' (p. 274) 

Lemma 13 Let a € Aiiji), i, ii, . . . , G {1, . . . , . Then aa,i defined as in 
il3\) satisfies 

< CNXo{x)-'^''+^^ < CNmax{M-'^''+^\l), (15) 
for all X in R^ . 

The lemma shows that the functions a^^i are in C^{R^). Together with 
(fT4|) this immediately implies that the vector fields i — 1, . . . , N have finite 
formal degree no greater than p. Just like identity (|T^ the following corollary 
is a trivial generalisation of Corollary 18 in [2], the result is also implicit in 
Kusuoka T (Proposition 14) 

Corollary 14 Suppose the vector fields (Vi, i — 0^ . . . ,d) satisfy the uniform 
Hormander condition. Then for any j > 1 there is a constant C2 > indepen- 
dent of f and t such that 

<C2t-0-lW2||V/||oo 

for allte (0,1], feC^{R^). 

We point out that the constant C2 does (via the constant M in the Hormander 
condition) depend on the underlying family of vector fields Vi . 



OX^^ . . . OXi^ 



sup 

ii,...,i,efl N\ 



d 



d 



■Ptf 
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4 Application to cubature on Wiener space 



4.1 The reduction operation 

In the iterated KLV method (Section [3]) the total error Ej) over the interval of 
approximation [0, T] is bounded by the sum of the individual errors Eg. over 
smaller time intervals. The KLV method is sequential. Starting with a unit 
mass particle at a single point in space time the measures evolve through time 
by replacing each particle at time ti with a family of particles at time U+i. 
Together these new particles have the same mass as their parent particle and 
are carefully positioned to provide a high order approximation to the diffusion 
of the underlying SDE. The algorithms introduced in section[5]can be used very 
effectively to perform a global redistribution of the mass on the particles alive 
at time ti so that an essentially minimal number of particles has positive mass. 
At the same time we do not increase the one step errors Eg- significantly or 
affect the order of the approximation. In this way we obtain (see Section 14.21) 
a global error bound over [0, T] for this algorithm that is of the same order (in 
the number of time steps) as the unmodified KLV method. On the other hand 
the blow up in the number of particles is radically reduced. 

The property of the intermediate measures we are targeting is to integrate 
Ptf correctly. To approximate the integral of a smooth function such as Ptf 
with respect to a discrete measure we need to find uniform functional approxi- 
mation schemes that apply to smooth functions on the support of this measure. 
By definition smooth functions can always be well approximated on balls by 
polynomials. However, only after one has set a fixed error bound e and a degree 
for the polynomials the size of the balls on which the approximation holds be- 
comes clear. The main idea will be to localise the intermediate particle measures 
Q into measures Qi, where each Qi has its support in such a good ball. We then 
replace^- using the algorithms of Section [5]- the measures Qi by reduced mea- 
sures Qi that integrate polynomial test functions of degree at most r correctly. 
In that way one knows that for a smooth function g 

J gdQi 

i 

is a good approximation to J gdQ.We subsequently prove that we can choose 
the localisation of the measure Q in a way that ensures that we increase the 
overall bound on the error of the approximation only by a constant factor and 
examine how well we can cover the support of the intermediate measures Q by 
balls for the localisation. 

A main idea for estimating s is to consider Taylor expansions of the function 
Ptf. We define p to be the minimal integer k such that the vector fields {Va, a € 
Ai{k)} uniformly span at each point of a; G (as in the UH condition). 
For g a smooth function on R^ let dg : R^ — > Hom(R'^ , R) denote the full 
derivative of g. The second order derivative d^g is then mapping 

i?^ ^ Hom(R'^, Hom(R'^, R)) = Hom(R^ ® R^, R). 
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The higher order derivatives can similarly be regarded as sections of 

Hom((RN)«\R). 



We define the r th degree Taylor approximation of g centred at a;o G to be 



Tayr{g,xo){y) = ^{d'g){xo)- 



(16) 



i=0 



and the remainder Rr{g,xo){y) by 

Rrig, xo){y) = g{y)- Tayr{g, xo){y). 

It is clear that the rth degree Taylor approximation centred at xq is a poly- 
nomial of degree at most r. Given u > and y G let B {y, u) denote the 
Euclidean ball of radius u > centred at y. Our estimate for the remainder of 
the polynomial approximation is the following: 

Lemma 15 Let t G (0,1]. T/ie remainder function Rr{Ptf ,XQ){y) is uniformly 
bounded on B (zq, u), i.e. 



\Rr{Ptf,Xo)\ 



< Ca 



,r+l 



iv/ll. 



where C4 — C2C3 is a constant independent of f , u and t. 
Proof. By Taylor's theorem we have ioi y E B {xq, u) 



\RriPtf,xo)iy)\< 



(r + 1)! 



\y - xoi 



r+1 



and we note that 



9\\^<C3ir,N) sup 

iiH hjiv=»"+l 



dxY 



Ptf{y) 



for some constant C3 that only depends on r and N. From Corollary [14] we see 
that 

nil Qijv 

-Ptf 



sup 

iiH |-ijv='"+l 



dx\' 



dx'^ 



< C2t-''P/^\\Wf\\ 



where C2 is the constant from Corollary [14] and the claim follows. ■ 

The bound on the remainder of the Taylor expansion of Ptf implies that 
cubature measures which integrate polynomials up to degree r correctly provide 
good approximations provided the support of the measure we are targeting is 
contained in a sufficiently small patch. 
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Proposition 16 Suppose the uniform Hormander condition is satisfied. Let 
t € (0,1] and ^ be a positive measure on with finite mass v satisfying 
supp{ii) C B{xo,u) for some u > 0, xq . Suppose a measure ft is a 

degree r cubature measure for fi (a reduced measure with respect to fj, and the 
polynomials of degree at most r) . Then 

\E,Ptf - E~,Ptf\ < CiV-^\\\7f\\^, 

where C4 is the constant from Lemma \J5\ and independent oft, /, x'o and u. 
Proof. We have 

E^^Ptf - E^Ptf - 

{E^ - E~^){Tayr{Ptf.xo)) + E^Rr{Ptf,xo) - Ef,R^{Ptf,xo) 

Since fi is a cubature measure and integrates polynomials of degree at most r 
correctly the first term of the sum vanishes. Lemma 1151 gives us the required 
bounds on the remaining terms. ■ 

Let /i be a discrete probability measure on R^ and {UjYj^^ be a collection 
of balls of radius u on i?^ that covers the support of /i. Then there exists a 
collection of positive measures fij, 1 < j < ^ such that /i^ _L /i^ for all i ^ j (i.e. 
the measures have disjoint support), 

e 

i=l 

and supp(^j) C Uj n supp (fi) . We call such a collection (^Uj,fij) a localisation 
of /i to the cover {UjY^^^ and say u is the radius of the localisation. 

Definition 17 We say that a measure Jl is a reduced measure with respect to 
the localisation (Uj, ^j'j and a finite set of integrable test functions P if there 

exists a localisation (Uj,]ljY_._^ of Jl such that for 1 < j < ( the measures Jlj 
are reduced measures ( see Definition with respect to /i^ and P. 

Note that the localisation of the reduced measure fi is with respect to the 
same cover as the original measure /i. It is trivial to show that reduced measures 
/X exist for any localisation (?/,■, /i ) . of a discrete probability measure and 
any finite set of integrable test functions P. Moreover the number of particles in 
the support of Jl is bounded above by (card (P) + 1) £. The following corollary is 
an immediate consequence of Proposition [TBI Let P in the following be a basis 
for the space of polynomials on R^ with degree at most r. 

Corollary 18 Let t < 1 , ji be a discrete probability measure on R^ and 
(Ujj^jj . a localisation of radius u. If ^ is a reduced measure with respect 
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to (Uj, ^j) and P we have 

\E,Ptf - Ef^Ptfl < Ci-^\\Vf\U 

where C4 is the constant from Lemma \T5\ and independent of t, f , u and the 
localisation of radius u. 

We define the Kusuoka-Lyons-Victoir transition KLV over a specified time 
interval [0,s], based on the cubature on Wiener space approach, and already 
used in the iterative method in Section [31 The transition KLV takes discrete 
measures on to discrete measure on and may be interpreted as a discrete 
Markov kernel. Given a measure = I'-i^xi on i?^ the new measure is 

obtained by solving differential equations along any path in the support of the 
cubature measure 

n 
i=l 

starting from any particle in the support of /i. We define 

I n 

KLV{p,s) = X!Xl'^J•^'^*=.-,(".)■ 
We are ready to consider recombination for the iterated KLV method. Let 
2? be a fc step partition = < < ... < tk = T of [0, T] the global 
time interval of the approximation and recall that sj = tj — ij-i. We also let 
u = {u2, . ■ . , Ufe-i) G R''~'^ where each uj > 0. Let P be a basis for the space of 
polynomials on R^ with degree at most r. For each time step Sj we first apply 
the KLV method to move particles forward in time to a measure Q. We then 
localise the measure Q and use the algorithm of Section [5] to compute a reduced 
measure with respect to the localised measure and replace Q by this reduced 
measure. The Uj determine the radius of the balls in the localisation of the 
measure in the j th iteration of the method. The polynomials in P serve as the 
test function in the reduction. 

More precisely we define two interrelated families Q^y^ (x) and (x) of 
measures. As base case we have the measures obtained by applying twice the 
KLV operation starting from the point mass at x. 

Q^^l (x) KLV{S., si) gg^„ (x) KLV{Q^^l (x) , s^) (17) 

For the recursion the measure (5p''„ {x) is defined to be a reduced measure with 
respect to any fixed localisation {UjiQj^u measure (x) with 

radius Uj and the set of test functions P (polynomials of degree at most r). We 
define Q^-p^^^ {x) by the relation 

Qg+i) (x) KLViQ^il (^) ' ^^+1) (l^) 
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for all z = 2, . . . , fc ~ 1. Note that we do not recombine after the first and 
last application of the KLV operation. The reduced measures Qp^„ (a;) are not 

unique even after we fix a localisation of Qp^„ (x) and a reduced measure may 
be computed using the reduction algorithms of Section [2j 



The main result of the section is the following: 

Theorem 19 For any choice of localisations (^Uj,Q^j^^{x)j^ with radius 

and any reduced measures Q^\{x) with respect to {Uj,Q^^^{x), 
functions P, 2 < i < k ~ 1, we have 



and test 



Ev.k ■= sup 



Prfix) ~ 



,0- + l)/2 



k — l m+1 

<\oat){si^^+j:j:j?^. 

i—l j—tn 



r+1 



CAT) 



fe-1 

^ (T - ij)'-p/2 



IV/llc 



(19) 



where Ci (T) and C5 (T) are constants independent of f and the choice localisa- 
tions with radius Ui.The constant C5 (T) can be taken equal to C4 ifT — ti < 1. 



Proof. The global error is bounded by 



Prfix) 



E (k) , N / 



< 



Prfix) 



fe-1 

E^w ,,^,PT-t J - E^i2-, (^j^V-ts/ kgO' {x)PT~tJ - E^ir, (^)Pr-tJ 



MY 



3=2 



'A(j) i^\Pr~tif - 



fc-1 

The first two terms and the terms in the second sum are the errors introduced 
by the KLV operation and can be bounded as in the proof of Theorem [9] 
The terms in the first sum may each be bounded by using Corollary 1181 ■ 

The bounds for the error derived in this section assume that the function / is 
Lipschitz. If / has more regularity it is clear different estimates can be applied 
to estimate the derivatives of Ptf giving alternate bounds for iJ^fe .Clearly a 
smaller number of balls in the localisations of the measures Q^^^ {x) reduces 
the computational complexity of the method. We have not discussed yet how 
to choose the localisation and the degree r in the reduction to optimise the 
computational complexity of the method (see Section |473)) . 
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4.2 Examples for the rate of convergence of the recombin- 
ing KLV method 

In this subsection we consider some particular choices of parameters for the 
recombining KLV method and examine their rate of convergence. We first fix 
for the remainder of this section a (family of) partitions T) for the time interval 
[0,T]. We recall a family of uneven partitions from Lyons, Victoir [T^ which 
has smaller time steps towards the end and is given by 

For 7 > TO — 1 the results in [T^] (see also Kusuoka [5]) show that 

5(j + l)/2 

(T - uy/^ 



k-lm+l U + l)/2 

T + T^^^TAln ^ ^« ^) (21) 

i—1 j—ni 



while for the case < 7 < to — 1 one obtains 

k~l m+l (j + l)/2 

(T - uy/ 



i—l j—ni 

In the following two examples we work with the partition defined in pop and 
the notation of Theorem [191 Using this particular choice of partitions ensures 
that the bound on the KLV error is of high order in the number of iterations 
k. 

Example 20 Lei 7 > to— 1, r — \m/p\ anduj — s^^^ °, where a := 2{[mjp ] +i] — 
0. Then 



sup 



fc-lm+l (j + l)/2 \ fe-1 Mni/p]p+l)/2 



< \^C, (T) 1^4/^ + (^e^ j + (^) E ^T-U)lmM.f2 j II V/ll 

< C8fc-(™-i)/'ri/'||v/|loo (22) 

where Cs = Cq (to, 7) {Ci (T) + C5 (T)) . 

Note that < p/2 — a < p/2 for all positive integers p and to and that for 
Sj < 1 we have Uj > Sj . In the next example we choose the radius of the balls 
in the reduction operation such that at each step in the iteration the bound on 
the recombination error matches the bound on the KLV error. 

Example 21 Let 7 > m — 1, m = r, i.e. the degree of the polynomials used in 
the reduction operation equals the degree of the cubature in the KLV method. Let 
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Uj, j — 2, . . . k — 1 be given by 



m+l 



Then 
sup 



m+l (j + l)/2 \ k-1 (m+l)/2 



<C9fc-('"-i)/2j,i/2||y^||^^ (23) 



w/iere Cg = Ce (m, 7) (Ci (T) + C5 (T)) , 



As before if T — ti < 1 the constants Cs and Cg can be taken to be 
Cq (to, 7) (Ci (1) + C4) .The parameters choosen in the above examples guar- 
antee high order convergence, but are not necessarily computationally optimal. 
In the following section we examine how for a fixed error e the choice of r 
and u can be varied to be closer to the optimal computational effort in the 
recombination operation. 

4.3 An optimisation 

This paper establishes stable higher order particle approximation methods where 
the computational effort involved grows polynomially with the number of time 
steps when the number of steps is large and the underlying system remains 
compact (see section |44|) . In concrete examples, an optimisation of the different 
aspects of this algorithm, under the constraint of fixed total error, leads to even 
more effective approaches - although we expect that different problems would 
benefit from different distributions of the computational effort. For example 
there is a trade off between the degree of the polynomials that are used as test 
functions and the size of the balls used to define the localisation of the measure 
for the recombination (smaller patches if we use higher degree polynomials in 
the test functions and we fix the error of the approximation) . 

Specifically, suppose we are given a discrete measure /i and the property 
we care about is the integral of against a smooth function g . As in our 
application to the KLV method we consider a reduced measure Jl (Definition [1]) 
with respect to the polynomials of degree at most r and a localisation of /x with 
radius at most S. The number of balls of radius 6 required to cover the support 
of jj, is at most of order (^) , where D is the diameter of supp(/j,) . Let e be 
the error of the approximation of J gdfi by J gdjl. 

Note that 

rr+l 

Cr+1 



e = 



(r + 1)! 
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for some c,.+i < En+...+««=r+i 
simple relation for S and r 



£(r + l)!V/('^+^) 



Cr+l 



Fixing the error e gives a 



(24) 



Let n be the number of particles in the support of /i. The computational com- 
plexity of the recombination operation as a function of 6, n and r is at most of 
order 



N 



r + N 
N 



\ogn 



r + N 
N 



which may be optimised subject to the constraint (|24p . 

Note that in our application to cubature on Wiener /x corresponds to (x) 
and the function g is given by Pr-tj /■ The calculation above also allows us to 
decide after each step of the iteration if it is of computational benefit to carry 
out a (full) recombination operation. 



4.4 Simple bounds on the number of test functions - cov- 
ering the support of the particle measures 

In this section we obtain upper bounds for the number of ODEs required to 
solve in the recombining KLV method with k iterations. For this it is sufficient 
to bound the number of balls in the cover of the localisations of the particle mea- 
sures uniformly for all k iterations. We first find a large ball B{x, p) that covers 
supp((5p^„ (x)), j = 1, . . . , fc — 1 and then estimate the number of balls that are 
required to cover B{x,p). The balls in the covers of the localisations will have 
to be sufficiently small to preserve the high order acuracy of the method. We 
can show that under the assumption that the vector fields Vi are bounded and 
satisfy the UH condition we have a high order method and the computational 
complexity is polynomial in k the number of iterations. Similar results can be 
obtained if the underlying system remains compact. 

The following theorem demonstrates that we can achieve the same rate of 
convergence in the number of iterations k as in Kusuoka's algorithm and the 
vanilla KLV method, but control the complexity of the method by an explicit 
polynomial in k. This compares to exponential growth in the vanilla KLV 
method without recombination, which despite its exponential growth leads to 
numerically highly effective algorithms (see e.g. Ninomiya, Victoir |13|). The 
estimates in this section are not designed to be optimal and can be improved. 
Closer to optimal choices for the radius Ui and degree r in the reduction op- 
eration have been discussed in Section 14.31 and may be used to decide if it is 
computationally efficient to recombine the particle measure at time U. 
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Theorem 22 Suppose the uniform Hormander condition is satisfied and the 
vector fields Vi are uniformly hounded by some constant M' > 0. We can 
achieve 



Ev,k = sup 



Prfix) - EQ^^>j,)f ^ Csk-^"'-'^/^T^/^Wf\U (25) 



while the number of test functions in the reduction operation, and hence the 
number of elementary ODEs to solve, grows polynomially in k. 

Proof. Let to > be the degree of the cubature in the KLV method. Fix 
the partition T) to ([201) fo'" some 7 > to — 1. As in Example let r = \m/p'\ 
and Uj = s^J"^ a = 2(rni/p]+i) 5^ the reduction operation. We note 



that the error E-D^k satisfies p5|) and it remains to show that the number of 
particles in the support of the measures Q^''^^ (a;) grows polynomially in fc, which 
is equivalent to the number of balls in the localisations growing polynomially in 
k. 

Note that if w e C^([0, l],R'^) is a continuous path of bounded variation of 
length L we have 

|a^-$i,xHI <M'L, 

where $ is the Ito functional defined in © , i.e. $1 is the point we obtain 
by solving the equation ([5]) along the path lo starting at x. Let L be given by 
L — maxi=i_....„length(ti;i) , the maximum of the lengths of the paths in the 
support of the degree to cubature formula on Wiener space over the unit time 
interval. Observe that by construction any particle in the support of {x) 
(compare the definition of the measures in (|18l) ) may be written as 

some ii,...,ij g {1, . . . , n}, the ujs,i are the rescaled paths defined in (0) and ® 
denotes to the concatenation of paths. For k sufficiently large we may assume 
Si <1 and we deduce that 

supp(g^^„ {x))<ZB ^x, M'L ^ s]''^ C B (x, M'LkT^''^ 

In the reduction operations we consider a basis of the polynomials of degree at 
most [m/p] and the measure is localised by balls of radius Uj which need to 
cover supp(Qp''^ i^))- For s.j < 1, i.e. for k sufficiently large, we have Uj > s^^^ 

and for our uneven family of partitions minj=2,...fc-i Sj^^ < s^^^ = (T /k'~'Y^'^ . 
Thus, the number of particles in each of the reduced measures is uniformly 
bounded above by times the number of balls of radius [T/k^Y^'^ 

required to cover the ball B [x, M' LkT^^-^) in N dimensional space, which is a 
polynomial of degree at most N {'yp/2 + 1) in k. ■ 

Similarly we can derive a result analogous to Theorem [22] if the underlying 
system remains compact. 
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